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Radiative energy transfer between closely spaced bodies is known to be significantly larger than 
that predicted by classical radiative transfer because of tunneling due to evanescent waves. The- 
oretical analysis of near-field radiative transfer is mainly restricted to radiative transfer between 
two half-spaces or spheres treated in the dipole approximation (very small sphere) or proximity 
force approximation (radius of sphere much greater than the gap). Sphere-sphere or sphere-plane 
configurations beyond the dipole approximation or proximity force approximation have not been at- 
tempted. In this work, the radiative energy transfer between two adjacent non-overlapping spheres 
of arbitrary diameters and gaps is analyzed numerically. For spheres of small diameter (compared to 
the wavelength), the results coincide with the dipole approximation. We see that the proximity force 
approximation is not valid for spheres with diameters much larger than the gap, even though this 
approximation is well established for calculating forces. From the numerical results, a regime map is 
constructed based on two non-dimensional length scales for the validity of different approximations. 



I. INTRODUCTION 



It is well known that thermal radiative transfer in the 
near-field is extremely different compared to classical 
radiative transfer in the far-field. Interference effects 
and, in particular, near-field effects due to tunneling of 
evanescent waves play a vital role. When the objects in- 
volved in energy transfer can support surface waves, the 
heat transfer can be enhanced by orders of magnitude 
compared to far-field values. The near-field exchange 
between two half-spaces has been well documented in 
literature P, H, H, 0, H, H, 0, Hi ■ Theoretical analysis of 
heat transfer between a sphere and a plane or between 
two spheres are limited to the sphere being approximated 
by a point dipole [1, d, [13, EH ■ In jg;] , the authors out- 
lined a method capable of dealing with near-field radia- 
tive transfer between a sphere and a fiat substrate when 
retardation effects can be neglected. However, numerical 
difficulties prevented them from obtaining a solution to 
the problem. It is possible to derive an asymptotic ex- 
pression for radiative transfer between two large spheres 
separated by a very small gap (gap is very small com- 
pared to the radius of either of the spheres) from the well 
known results of radiative transfer between two semi- 
infinite objects This idea is used extensively in de- 
termining van der Waals or Casimir force between macro- 
scopic curved objects and is known as the proximity force 
approximation [Tl, [l^.[T5l [lB| and expressed as the prox- 
imity force theorem [l7l|. The usage of the proximity force 
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approximation to determine near-field radiative transfer 
between curved surfaces has not been verified by other 
numerical solutions or by experiments. Generally, the di- 
ameter of the sphere involved in the experiments range 
from a few microns to a few tens of microns, which is 
no longer in the point dipole approximation. Hence a 
complete scattering solution to the problem is necessary 
to verify with experiments as well as gauging the va- 
lidity of simpler models. In this paper we investigate, 
for the first time, the near-field radiative heat transfer 
between two spheres using the dyadic Green's function 
(DGF) of the vector Helmholtz equation [l 8|. [ l9l. [20| and 
the fluctuation-dissipation theorem [2lll22l.l23j. This for- 
malism of fluctuational electrodynamics was pioneered by 
Rytov [i^] and is used widely for analyzing near-field ra- 
diative transfer 

S II i, i, 0, Hi . Though the emphasis 
of this work is on near-field radiative transfer, it should 
be pointed that this formalism includes far-field contri- 
bution to radiative transfer. 

Electromagnetic scattering by a sphere has been very 
well studied since the seminal work of Mie, almost a cen- 
tury ago. The two sphere scalar and vector scattering 
problems have also been investigated for almost the same 
amount of time by many authors [1^, [l^, [l^l . The two 
sphere problem involves expansion of the field in terms of 
the vector spherical waves of each of the spheres and re- 
expansion of the vector spherical waves of one sphere in 
terms of the vector spherical waves of the second sphere 
in order to satisfy the boundary conditions. The two 
sphere problem, and multiple sphere scattering in gen- 
eral, is especially tougher due to the computational de- 
mands of determining translation coefficients for the (vec- 
tor) spherical wave functions [H, [28|, |29|, [s^l . Recurrence 
relations for the scalar [3l| and vector spherical waves [s^l 
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have reduced the computational complexity considerably. 
In this work, the DGF for the two sphere configuration 
is determined by satisfying the boundary conditions for 
fields on the surface of the two spheres. The translation 
coefficients are determined using the recurrence relations 
in [3l|,[32. 

The paper is arranged as follows. In Section |TT1 sim- 
plified, asymptotic results for the radiative heat trans- 
fer between two spheres of equal radii, based on the 
dipole approximation and proximity force approxima- 
tion are presented. In Section IIIII the DGF formulation 
and fluctuation-dissipation theorem are introduced. The 
Poynting vector is expressed in terms of the DGF and 
material properties. In Section IIVI the two sphere prob- 
lem is described and the DGF for this configuration is 
determined in terms of the vector spherical waves of the 
two spheres. In Section |Vl the expression for radiative 
flux, and thus the spectral conductance, from one sphere 
to another is determined. Details regarding the conver- 
gence of the series solution and numerical solutions for 
sphere sizes up to 20 fim in diameter is presented in Sec- 
tion IVll 



II. ASYMPTOTIC RESULTS FOR 
NEAR-FIELD THERMAL RADIATION 

The purpose of this section is to present an asymp- 
totic results for the radiative heat transfer between two 
spheres in the dipole limit as well as the proximity force 
approximation limit. It is possible to define a radiative 
conductance between the two spheres at temperatures Ta 
and Tb as: 



G= lim 

Ta^Ts \Ta-Tb\ 



(1) 



where G (units WK~^) is the radiative conductance, 
P{Ta,Tb) is the rate of heat transfer between the two 
spheres at Ta and Tb- It should be noted that G is a 
function of temperature. If the radii of the two spheres 
involved in near-field radiative transfer are much smaller 
than the thermal wavelength (At = hc/kBT w 7.63/im 
at 300 K) , the spheres can be treated as point dipoles and 
the conductance between the spheres is given by : 



G{T) 



dO (uj, T) „ , , , - , 

' a'i{io)a'i{u;)dw (2) 



where 8(w,T) = hw / {exp{huj /kBT) — 1); a'( {to) and 
Q!2 (to) are the imaginary parts of the polarizability of 
the spheres; d is the center-to-center distance between 
the spheres. This d^^ behavior of conductance between 
two spheres is valid only when d ^ Xt and d 3> i?i + -R2 
i. 

Calculating the forces or heat transfer between two 
curved objects as the size of the objects increases become 
computationally difficult. For spheres of large diameters 
and much smaller gaps, the proximity force approxima- 
tion is very useful in calculating the forces between the 



objects and used widely to calculate van der Waals or 
Casimir forces [H, [H, [H, [l3 with the knowledge of 
the same forces between two half-spaces. The proxim- 
ity force approximation was used to calculate the con- 
ductance between a sphere and a flat surface ^12] from 
the results of the radiative heat transfer between two 
half-spaces 0, 0, [H, [sj . The conductance per unit area 
for flat plates is the radiative heat transfer coefficient, 
h (units Wm~^K~^). The heat transfer coefficient h{x) 
in the near-field, especially when dominated by surface 
polaritons, is known to vary as l/x^, where x is the gap 
between the half-spaces [ll|, [11]. The spheres are sep- 
arated by a minimum gap x. We will assume that the 
spheres are of equal radii, R. The conductance between 
two spheres is computed by approximating the spheres 
to be flat surfaces of varying gap. By doing so we get a 
relation between G and h given by 



G{x;R) = TTxRh{x) 



(3) 



Since h{x) varies as it is expected that G{x;R) 

varies as l/x. A 1/x variation of conductance as 
determined by a more rigorous theory can be taken as 
evidence of the validity of the proximity type approxi- 
mation at the value of the gap. Despite their origin in 
fluctuations of the electromagnetic field, a significant 
difference between force and fiux is that the force 
decays to a negligible quantity in the far-field whereas 
thermal flux attains a finite value in the far-field. This 
implies that those parts of the spheres with larger 
gaps contribute very little to force whereas they could 
contribute significanly to flux because of the larger areas 
involved. Hence a proximity force type approximation 
would be valid only when the heat transfer is dominated 
by contributions from the near-field region. Therefore, 
we expect the result of Eq. [3] to be valid only when R 
is small enough that near-field radiation dominates and 
x/R 0. The discussion in Section IVTl shows that this 
is indeed true. 



III. ELECTROMAGNETIC FORMULATION 

All the materials are assumed to be non-magnetic 
and defined by a complex, frequency dependent dielec- 
tric function, e(w). To compute the radiative transfer we 
follow the method pioneered by Rytov [l^, |2^ in which 
the source for radiation is the thermal fluctuations of 
charges. The Fourier component of the fluctuating elec- 
tric field, E{ri,uj), and magnetic field, H{ri,uj), at any 
point, 7*1, outside a volume containing the sources is given 

by mm 

E{ri,uj) = iujfio d^rGeiri,r,uj) ■ J{r,Lu) (4) 
Jv 

H{ri,uj) = / dV^,,(ri,r,c^)- J(r,w), (5) 
Jv 
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where Ge{ri,r,uj) and Gh{ri,r,uj), the dyadic Green's 
functions due to a point source at r, are related by 

Gh{ri,r,uj) — Vi x Ge(ri,r,w); J{r,uj) is the Fourier 
component of the current due to thermal fluctuations; 
and fj,o is the permeability of vacuum. The integration 
is performed over the entire volume V containing the 
source. The DGFs themselves obey the following equa- 
tions 

V X V X We{r,r') - %(r)^e(r, r') = T<5(r - r')(6) 

where I is the identity dyad and 6{r — r') is the Dirac- 
delta function. At the boundary between two dielectric 
(possibly lossy) materials, the DGF satisfies the following 
boundary conditions to ensure continuity of tangential 

I 



The fluctuation-dissipation theorem states that the 
cross spectral density of different components of a fluc- 
tuating current source in equilibrium at a temperature T 
is given by [2l[: 

{Ji{r,uj)J„^{r = Si,n5{r~r ), 

TT 

(10) 

I 



where the superscript T stands for the transpose of the 
dyad. Once the Green's function for the given config- 
uration is determined, the above integral is computed 
numerically. Determining the DGF is not a trivial task 
and the next section is devoted to determining the DGF 
in the case of the two sphere configuration. 



IV. TWO SPHERE PROBLEM 

The configuration of the two spheres is shown in Fig. 
[TJ At the center of each sphere is a coordinate system. 
Without loss of generality, the two spheres are arranged 
such that the z-axes of both coordinate systems pass 
through the line joining the centers. The x-axis (and 
y-axis) of both systems are parallel to each other so that 
a given point in space has the same i/i coordinate value 



electric and magnetic fields: 

n X Ge(ri,r') = n X Ge(r2,r') (7) 

n X V X Ge(ri, r') = n X V X Ge(r2, r') (8) 

where ri and r2 are points on either side of the bound- 
ary and ri is a unit normal to the boundary surface at 
ri (or r2). In order to compute the spectral Poynting 
vector at ri, we must compute the cross spectral density 
of Ei{ri,t) and Hj{ri,t), (Eii^H*^), where the * denotes 
the complex conjugate, the brackets denote a statistical 
ensemble average, and i and j refer to the three carte- 
sian components (i ^ j). From Eq. we can write an 
expression for {Eii^H*^) as: 



(9) 

I 

where e"(w) is the imaginary part of the dielectric func- 
tion of the source, €o is the permittivity of vacuum, and 
0(w, T) is given by hiu/ {exp{hLj /ksT) — 1), where 27rfi, is 
Planck's constant and fc^ is Boltzmann's constant. Using 
Eq. dni) and Eq. (dH]), we have 



(11) 

I 

in both systems. The two spheres are at temperatures 
Ta and Tb- In order to determine the radiative transfer, 
we have to determine the DGF when the source point is 
in the interior of one of the spheres. We shall take the 
Dirac-delta source point to be in the interior of sphere A. 
The most convenient way of dealing with DGF in spher- 
ical coordinates is to expand the DGF in terms of vector 
spherical waves jis!] , which are solutions of 

V X V X P(r) - fc2p(r) = (12) 

The vector spherical waves we will need are given by [isj : 

M^^fcr) = z\^\kr)vl-^\e,cp) (13) 

'-^VWTT)v^\e,^) (14) 



(E,{ri,u)H;{r,,Lo)) = ^u;^lo d?r dV'{Ge„ (n, r, uj)Gl^^ (n, r', u){Ji{r, uo)r^{r' , lo))} 



{E,{n,uj)H;{ruu)) 



ieoe"iu)t^oto^e{io,T) 



d^r{{Ge{ri,r,uj).Gl*{ri,r,uj)),j} 
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FIG. 1: Two sphere configuration. Two non-overlapping 
spheres of radii a and b are separated by a distance d. At 
the center of each sphere is a spherical coordinate system ori- 
ented such that the two spheres lie along the common z-axis. 
The x-axes and y-axes are also oriented such that for a given 
location in space, the coordinate is the same in both co- 
ordinate systems. In this figure, the point P has coordinates 
(rA,6A,<t>A) and {rB,OB,(t>B) such that (jiA = 0s. Region 
A(B) refers to the interior of sphere A(B). Region C is the 
exterior of both spheres and is taken to be vacuum. 



where the integration domain refers to the surface of a 
sphere of unit radius and is a differential area element 
on such a sphere. The vector spherical waves Mi^(kr) 
and (kr) are related by N^^ (kr) = ^ V x Mj;^J (kr) 

and M^^ikr) = -^V x Nl^{kr). Any solution to Eq. 
[12] can be expressed as a linear combination of the vector 
spherical waves. Further properties of spherical harmon- 
ics and vector spherical harmonics useful for this analysis 
is included in Appendix [X] and Appendix A of [1^ . In 
this particular case, the field is a linear combination of 
vector spherical waves of the two coordinate systems as 
shown in Fig. [TJ To satisfy the boundary conditions on 
the surface of each sphere the vector spherical waves of 
one coordinate system should be expressed in terms of 
the vector spherical waves of the other coordinate sys- 
tem. This is what is achieved by means of translation 
addition theorems for vector spherical waves. 



A. Coefficients for translation addition theorems 

The vector translation addition theorem [3, US) 113 
states that: 



where M.l^{kr) and N\^l{kr) are vector spherical waves 
of order {l,m). I can take integer values from to oo. 
For each |m| < I. The superscript p refers to the ra- 
dial behavior of the waves. For p = 1, the M and N 
waves are regular waves and remain finite at the origin 
and zl^\kr) is the spherical bessel function of order /. 
For p = 3, the M and N waves are outgoing spheri- 
cal waves that are singular at the origin and z\ {kr) is 
the spherical hankel function of the first kind of order 



The radial function C/''^ (x) = (^xzfix) 



Ml 



■V, 



Wi 



Ira 



(9, (/)), and VJ^^(6', (/>) are vector spherical harmonics 
of order {l,m). The three vector spherical harmonics 
can be expressed in terms of the spherical harmonics, 
Yj„: (9, 6) as: 



.(3), 



V, 



(1) 



Irn 



Ira 



(15a) 



(2)/ 



V, 



(3), 



do sinO 



(15b) 



O^T^r- + (t}-—^Yi„i ] (15c) 
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sinO 



The vector spherical harmonics are orthonormal to each 
other and satisfy the following relation: 



mq 



(16) 



Ml'Jikr,)^ ^ [4-(^-fcd)M(«)(fcr„)^ 

I.L=-U 

iy = l 

Bi^^(+kd)Nif{kr,)\ 



^/™(fcn)= E [Bl^,{+kd)Ml''^{kray 

u = l 

Ai,™(+fcd)iV(^)(fcr,)] 



(17a) 



(17b) 



The position vectors and refer to the same loca- 
tion in space in coordinate systems A and B respectively. 
Computing the coefficients A^l^{+kd) and B\™{+kd) has 
been the topic of many publications 28, 29, 30, 32]. Gen- 
erally, the expressions for the coefficients require calcu- 
lations of Wigner 3j symbols which involve calculations 
of large number of factorials, making it computationally 
expensive. Recurrence relations for computing the coef- 
ficients efficiently have been proposed by Chew j32j. In 
the case of the two sphere problem, with translation along 
the z-axis alone, Eq. [T71 simplifies so that the coefficients 
are non-zero for fi = m alone. 

In the region C (exterior to both spheres), the elec- 
tric and magnetic fields should be expanded in terms 
of outgoing vector spherical waves of both coordinate 
systems so that waves decay as 1/r as r oo. Hence 



Mif^(kri,) and Nl^l{krb) need to be expressed in terms 

of Mi^{kra) and N^m{kra) on the surface of sphere A 
and vice versa. Since \ra\ ^ a < d for all points on 
the surface of sphere A, only the regular vector spheri- 



r(3) 



cal waves or Mlm{kra) and Nlm{kra) should be used 
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In addition to A^J^{+kd) and B^™ (+fc(i), we will also 
need A^^{—kd) and i^i™ (—kd) which can be obtained 
through symmetry relations [38|. For further details re- 
garding the computation of the recurrence relations, the 
reader is referred to [3l|, [12] . 



infinite medium, Go and one that takes into account the 

scattering, Gsc- In this case, the source point is confined 
to the interior of sphere A. The DGF for source point 
in sphere A, assuming the whole space to of the same 
material, is given by: 



B. DGF - vector spherical wave expansion 

The DGF for any configuration can be split into two 
parts - one that corresponds to a Dirac-delta source in an 



l = oo 
m = l 

E 

1 = 1 



M^^ik,ra)M!XikaK) 



N^^ikara)N^:iJk,K 



?i) 



AkaK) 



if ra 
if ra 



(18) 



In particular, we are interested in the case where ra > r'^ 
since the source is inside the sphere A whereas the bound- 
ary of interest is the surface of the sphere. The part of 
the DGF that depends on the boundaries takes differ- 
ent forms in the three regions. A, B, and C. Inside A, 
the DGF is a combination of Go and Gsc, whereas out- 
side A the DGF is entirely Gsc- Each term, Mi^{kara) 



or Ni^{kara), in Eq. [TH] can be thought of as an in- 
dependent vector spherical waves that produces scat- 
tered waves, i.e. coefficients of scattered waves due to 
^hnikaTa) ipv N^^ (kaVa)) are Completely decoupled 
from vector spherical waves of other orders. Let us con- 
sider the scattered field due to {kaVa). The scat- 
tered field in the three regions is given by: 



oo 

E 

v'—[rn,l) 



ika[{AiZM^]l{kara) + A^^MJ^ (fcar,)) + (B^^M*^ (fc.rfc) + B'ZNi'rIikan))] , in A 
ikf[{CZMl^Aikfra) + CiZNi^likfra)) + {Dif^Mi''l{kfn) + Di^^Ni'llikfn))] , in C (19) 
ih[{EZMi,^2ihra) + E'^Ni^lihra)) + (i^^^M^ (fc^rt) + F'ZNi'^ihn))] , in B 



where the symbol (to, 1) refers to the greater of m and 1. 

is a coefficient of a vector spherical waves of order 
[v, m) that it is produced by a vector spherical waves of 
order (Z, to) . The superscript M (N) is to indicate that 
it is a coefficient of a M (N) wave. In practice, the 
upper limit for the summation is resticted to a value Nm 
which depends on kfd. The appropriate value of Nm will 
be discussed in Section |Vll Using Eq. [191 Eq. [3 and Eq. 
m the following set of coupled linear equations can be 
obtained for the coefficients of the vector spherical waves 
in the scattered field in region C: 



clZ + M<^) E Ks-(-fc/d)-H 

i/=(ma) (20b) 

DZA^;;;,{-kfd)]=Q 



N„ 



CZB;Zi+kfd)] =0 



clX + M'^) E [DZA;'^i-kfd)+ 

v—{'m^l) 

DZ.B'^^Zi-kfd)] ^p^5,i 



(20a) 



KiB'^zi+kfdy 



" \J^n^,' ^"^ ""^ ^ ■ (20d) 

cZA;^{+kfd)] = 

where rj ranges from (to, 1) to Nm- Uriia) and Vrj{a) are 
Mie coefficients that one encounters in the scattering of 
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a spherical wave by a single sphere and are given by 

^ _ fcaC^^^(fcaQ)4^Hfc/a) - kfC}^\kfa)zl^\kaa) 

kaC}q\kaa)z''i^\kfa) - kfCip{kfa)z^rl\kaa) 

(21a) 

ka(!;^'^\kfa)zl^\kaa) ~ kfC!;^'^\kaa)zl^^\kfa) 

kaC,'i'\kfa)zl^\kaa) - kfC^^\kaa)zi^^\kfa) 

(21b) 

Expressions for (b) and {b) are obtained by replacing 
ka and a by fcb and b respectively. If the original wave 

is ^imi^afa) instead of Ml^{kara), the only difference 
is that the right hand side (RHS) of Eq. I20al becomes 
and the RHS of Eq. I20bl becomes 5ni ■ _p*^ and are 
given by 



-i/{kfa) 



kaaC}q\kaa)z'i^\kfa) - kfaC^^\kfa)z.^^\kaa) 



(22a) 



i/{kfa) 



kaaC^ri'\kfa)z'^l\kaa) - kfaC}^\kaa)z':i^\kfa) 



(22b) 



For a given m, we have {N„i — (m, 1) + 1) M waves and 
{N„i — (m, 1) + 1) AT waves. We see that the left hand 
side (LHS) for a given value of m remains the same while 
the only difference is in the RHS (p^^ and ) . Once all 
the coefficients in Eq. [50] are obtained, the DGF due to 
scattering, Gsc, and its curl, V x Gsa can be written as: 



Z,L'=(1,T^) 



{DZM!.'J,{kfn)- 

r(3)/ 



-^'IN 7\r(3) / 



{CZMliiikfr,) + CZNi.-:>{kfr,)) + 

{kfn)) 

{C'J^M!,tl{kfra) + C'J^N^tl{kfra)) + 
{D'l^;^Mlt(kfn) + D'ZNi?l{kfn)) 



N, 



(1) 



(kaK) 



(23) 



V X Gsc{ra,K) = ik] 



E 

m=~N„ 



i-iy 



( ^Im ^I'm {kfTa) 



{Di^N^J'likfn) 4 
/^iV(^(fc/r<,)4 



IN 



(3) 



(kfTa))- 



(3), 
'm \ 

(3), 



i/m(fc/r6)) 



C'J^Ml^ilikfra))- 



ivii/rn 



{kfn)) 



Mi'LikaK) 



(24) 



V. RADIATIVE FLUX 

The radiative heat transfer between the two spheres is 
calculated from the Poynting vector normal to the surface 
of sphere B, which in turn depends on the tangential 
fields on the surface of sphere B. The expression for the 



DGF can be modified to reflect tangential and normal 
fields on the surface of sphere B by using Eq. [20c] and Eq. 
[20dlin Eq. [23] and Eq. [2l and eliminating M;^^(/c/ra) 
and N^mikfra) to result in the following equations: 



Gsc{ra,K) 



b ^ 

i.,. = (l.m) 



(-ir 



Dl't^Y\hi.r„) 

D'SL(,\'\k^r,) 

ktfic}^-'' (kjh)z'-J-'\k^h)-kjhc}J-\k^h)z'-J-'^ (kjh) 



Vii^^iOhAh) + 

v^l\{ehAh) + 



(25) 



(fc,/fcj)D^^zp'(fcbrfc)(^,y(,.+l)/fc,b) 

+ [similar terms with primed coefficients] N^^^^^kar'^) 
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m = Nm 



Z,z. = (l,m) 
m.= -N.,-n 



ik,/kf)Dll4^Hk,r,) 

ktbcl'^\kfb)zl^''(khb)-kfb<;l'^\k,,b)zl^\kfb) ^ 

Dll4^\{kbrb)iy/i'{i'+l)/ktb) 
kbbcl^\kbb)zl^\kfb)-kfb<;l^\kjb)z''J-\kbb) 



V}?rli0b,M + 



[similar terms with primed coefficients] _^^{kar'^) 
I 



(26) 



Using Eq. [TTl Eq. [^H Eq. [211 Eq. IA2i Eq. IA3[ Eq. IA4[ expression for the spectral radiative transfer between the 
Eq. IA5[ and some algebraic manipulation to yield this two spheres, one at temperature and the other at Tg: 



P{u;;Ta,Tb) = (e(w,TA) - e(w,ri3))-x 



m,l,f3 



x^jib) 



zf\kib) 



r['\k^a)D'^^ 



4^ikfb) 



zl'-\k^a}D'/'„ 



rfHkfb) 



ri'\kfa)DX'i 



rjs'kkfb) 



(^) 1-^(^)1'- 



(27) 



r 



where (a) 
kfaC,l^\kfa)zl^\kaa) 



kaaC!;j^\kaa)zrl\kfa) 



and 



kaaC^^\kfa)zil\kaa) ~ kfa(^,j^\kaa)zlj^\kfa). 



Just 



as in Eq. [U it is possible to define a spectral radiative 
conductance between the two spheres at a temperature 
T {Ta^T,Tb^T) as: 



G (w; T) = lim 



P{uj-Ta,Tb) 



E 



Ta^Tu \Ta-Tb\ 

(^) 



= kf 



2„X 



zl^\k^a)D'^''l 



x^{b) 



zjpikfb) 
r[^\k^a)D'l 

■/3 



Z^^\kfb) 



(e^ - 1)' b 

2 



z[^\k,a)D\^„ 



rf\kfb) 
rl'\kfa)D'i,'' 



(^) 1-^(^)1' 



+ 



(28) 



where X = huj/kBT. It can be seen from Eq. \7E\ that 
the spectral conductance has units of kB{JK^^) and can 
split into two parts, one that depends on temperature 
and the other that is obtained from the DGF of the two 
sphere problem. The radiative conductance between the 
two particles that one would measure in an experiment 



is the integral of G {uj;T). 



Gt{T) = J G{Lu;T)dLo 





(29) 



8 



VI. NUMERICAL RESULTS 

To ensure that the program written to determine the 
near-field radiative transfer heat transfer is not misbe- 
having, a few tests can be performed. One of them is 
agreement between the numerical results and the ana- 
lytic expression for conductance in the point dipole limit. 
For spherical particles in the 1 nm to 50 nm radius, the 
numerical results agree well with the expression for con- 
ductivity in the point dipole limit [33 | . In addition to 
this another test to ensure the correctness of the method 
is based on the principle of detailed balance. The radia- 
tive conductance between two spheres of arbitrary radii 
must be independent of the numbering scheme for nam- 
ing the particles, i.e G12 — G21, where the subscripts 
1 and 2 refer to the two spheres. This is necessary to 
ensure that when the two spheres are at the same tem- 
perature the net heat transfer between the two spheres is 
zero. It is indeed seen from results shown in Table [J that 
by switching the position if the spheres, keeping the gap 
the same, results in the same value of conductance (the 
relative errors are generally of the order 10^^''') 

TABLE I: Conductance obtained for spheres of unequal radii. 
By swapping the radii of the spheres, it is seen that the value 
of conductance remains the same. 



Gap (/xm) 


a {fJ.m) 


b (pm) 


Conductance (WK"^) 


0.5 


1 


2 


1.63816 xlO"" 


0.5 


2 


1 


1.63816 xlO"" 


0.8 


2 


3 


3.34409 xlO"" 


0.8 


3 


2 


3.34409 xlO"" 



A. Convergence analysis 

Though the two-sphere scattering problem has been 
discussed in literature, the near-field interaction between 
the two spheres has not been analyzed in detail. In par- 
ticular, the number of terms required for convergence for 
the near-field problem has not been mentioned at all. 
For Mie scattering by a single sphere of radius a and 
wavevector magnitude k the number of terms for conver- 
gence, Nconv, is given by [H: 

Nconv = l + ka + 3(/ca)i/3 (30) 

For the two sphere problem a slightly different criterion 
based on the center to center distance between the two 
spheres is proposed and given by [40l |: 

Nconv = ^ekd (31) 

where e is the base of the natural logarithm. Both Eq. 
1301 and Eq. [3T] are valid criteria for computing far-field 



qualities, like the scattering coefficient. For the near- 
field problem, we expect the gap between the two spheres 
to be of great importance. To determine the number of 
terms required for convergence of near-field quantities we 
seek parallels to the much simpler and often investigated 
problem of near-field transfer between two half-spaces. 
In the near-field two half-space problem, the equivalent 
of the number of terms for convergence is the truncation 
for the in-plane wave vector. For a gap x between the 
two half-spaces, the predominant contribution to radia- 
tive flux in those frequency regions where electromagnetic 
surface waves are important is from in-plane wave vectors 
up to the order of kin ~ l/s. This is true when \/|e| ~ 1, 
where e is the dielectric function of the material of the 
sphere. This dielectric function of silica between the fre- 
quencies of 0.04 eV to 0.16 eV satisfies this condition, in 
which range the maximum value of -y/jij is « 2.1. If we 
can draw an analogy between the in-plane wave vectors 
in the two half-space problem and the two sphere prob- 
lem, we can propose a criterion for convergence. In the 
case of the two-sphere problem, the in-plane wave vector 
equivalent is given by the wavelength of periodic varia- 
tions on the surface of the sphere. A given vector spheri- 
cal harmonic, Vil^{9, (f), (j> = 1, 2, or 3) corresponds to a 
variation exp(im(/)) along the equator of the sphere. The 
period of this particular vector spherical harmonic (along 
the equator) is 2'Ka/m and the corresponding wave vector 
is m/a, where a is the radius of the sphere. For a given 
gap X between two spheres of radii a (equal for now) , the 
number of terms Nconv for convergence should be cho- 
sen such that {Nconv/a)x ~ 1. Hence the convergence 
criterion for near-field effects is given by: 

a 

Nconv ~ (^^) 
X 

For spheres of unequal radii the number of terms depends 
on the greater of the two radii. The convergence crite- 
rion proposed here is an upper limit and depending on 
the optical properties of the spheres, it could be consid- 
erably lesser. For spheres which exhibit surface wave res- 
onances, Eq. [32] is a valid convergence criterion as shown 
by our numerical investigations. Depending on the con- 
figuration of the spheres, the criterion for convergence is 
given by max{^ekd, ^). 

Since silica spheres are easily available for experimen- 
tal investigation, we will present our numerical results 
of radiative transfer between two silica spheres of equal 
radii at 300 K. In addition, silica is a polar material and 
hence supports surface phonon polaritons in certain fre- 
quency ranges. The dielectric function of silica is taken 
from [41[. The real part of the dielectric function is neg- 
ative for silica in two frequency ranges in the IR - from 
0.055 to 0.07 eV and 0.14 to 0.16 eV. It is expected that 
(and shown later) that surface phonon polariton reso- 
nances occur in the these frequency ranges. To confirm 
the prediction of Eq. [321 the contribution to the spectral 
conductance from each (/, m) mode is plotted as a func- 
tion of I for different values of m in Fig. [21 The spectral 
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conductance is plotted at two frequencies - 0.061 eV and 
0.045 eV. The surface polariton mode exists at 0.061 eV 
but not at 0.045 eV. 



2.0 10 



m = 0, 0.045 eV 

m = 1, 0.045 eV 

m = 1 , 0.061 eV 
/PI = 0, 0.061 eV 




a = b = 1 )im 
d = 20.1 )im 



§. 0.0 10° 



> 1.2 10" 



.9.0 10 



10 



o 

I 6.0 10-^° 

T3 

^3.0 10-'° 
o 

(§■ 0.0 10° 



a = b = 1 )im 
frequency = 0.061 eV 
kd = 6.21 



gap = 100 nm 
y gap = 200 nm 



50 100 150 200 250 
/ index 



FIG. 3: Plot of spectral conductance at 0.061 eV between two 
spheres with equal radii of 10 /im at gaps of 100 nm and 200 
nm. The curves plotted are for m = 0. 



FIG. 2: Plotted in this figure is the contribution to spectral 
conductance as a function of I for m = 0, 1 at 0.061 eV and 
0.045 eV for two spheres of equal radii a = 6 = 10 fim at 
a gap of 100 nm (d = 20.1 /im). Surface phonon polaritons 
contribute significantly to the radiative transfer at 0.061 eV 
and not at 0.045 eV. Hence the number of terms required for 
convergence is significantly lesser than that prescribed by Eq. 



In Fig. [3] the spectral conductance at 0.061 eV for m 
= between two spheres of radius 10 is plotted as 
a function of I for two different gaps - 100 nm and 200 
nm. As the gap doubles, the number of terms required 
for convergence approximately reduces by half, confirm- 
ing Eq. [321 For instance, the spectral conductance is 
1.02 X 10-"WK"^eV"^ at I = 181 for 100 nm gap and 
I = 99 for 200 nm gap. The contribution to spectral con- 
ductance from smaller values oil {I ^ 15), corresponding 
to propagating waves, does not change much with vari- 
ation in gap. It is the contribution from larger values 
of I, corresponding to surface waves, that changes appre- 
ciably with gap. Based on these results we use at least 
Nmax = '2a/{d — 2a) terms in our computations. Be- 
cause of the computation difficulties in solving Eq. [20l 
we present results for a maximum radius to gap ratio of 
100 for computations at one frequency and 100 frequency 
points (unequally spaced, with greater density in the res- 
onant parts of the spectrum) over the spectrum from 0.04 
eV to 0.16 eV. Equation [201 is solved using the software 
package Mathematica. 

As mentioned earlier, the symmetry associated with 
the two-sphere problem allows for solving for the contri- 
bution from each value of the m independently, starting 
from TO = and proceeding with increasing values of 
TO. As TO increases, the contribution to conductance de- 
creases as shown in Fig. [H The computation is stopped 
when a vacule of to is reached such that the contribution 
to conductance is less than 5 x 10"'^ times the contri- 
bution from TO = 0. Even though the contribution to 



conductance is significant for terms with I « Nmax, the 
contributions from to drops much faster. This is for- 
tunate - the time taken to compute the conductance is 
proportional to the number of values of to required. 

2.5 



3. 2.0 

CD 
O 

I 1.5 

T3 

o 1.0 
o 

^ 0.5 

CD 
Q. 
C/3 



0.0 



d = 200 nm 



d = 400 nm 
a = 20 )am 








5 10 15 
m index 



FIG. 4: Contribution to spectral conductance from each value 
of m at 0.061 eV. The curves shown are for spheres of radius 
20 /im and gaps of 200 nm and 400 nm between them. 



B. Spectral conductance 

Unlike the case of near-field radiative transfer between 
two half-spaces, where the conductance is a function of 
only the gap between the half-spaces (and the optical 
properties of the half-spaces and interveing medium), 
the conductance in the case of sphere-to-sphere radia- 
tive transfer varies as a function of the gap as well as the 
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size of the sphere. A gradual transition occurs from a re- 
gion of near-field dominated radiation to that of far-field 
dominated radiation. In Fig. [5] the spectral conductance 
between two silica spheres of 1 /im radius is plotted as a 
function of frequency. We see from the two peaks that 
the conductance is dominated by the surface phonon po- 
lariton regions. 




d = 50 nm 



d = 200 nm 



d = 100 nm 



0.08 0.12 
Frequency (eV) 




0.16 



FIG. 5: Plot of spectral conductance between two silica 
spheres of 1 radii at gaps of 50 nm, 100 nm, and 200 
nm from 0.04 eV to 0.16 eV. Surface phonon polaritons in 
the 0.055 to 0.07 eV range and in the 0.14 to 0.16 eV range 
contribute to the two peaks seen in the figure. 



The spectral conductance between the two spheres for 
larger radii, shown in Fig. [6l displays several features 
of interest. The ratio of radius to gap is maintained the 
same as in Fig. [51 Though the height of the peaks remain 
approximately the same in this figure as well as Fig. O 
the significant difference is from the contribution from 
those frequency regions that do not support surface po- 
laritons (0.04-0.055 eV, 0.07-0.14 eV). The contribution 
to the conductance from these ranges do not vary signif- 
icantly with gap (as long as gap <^ radius) . The relation 
between the ratio of gap to radius and the contribution 
to spectral conductance is also illustrated in Fig. [71 In 
Fig. [71 the spectral conductance of spheres of radii 1 
/im, 2 /im, and 5 at gaps of 100 nm, 200 nm, and 
500 nm respectively are shown. Since the ratio of gap 
to radius is a constant (0.1), we expect from the asymp- 
totic theory that value of spectral conductance should 
also be the same in all three cases. We see from the 
data that this is approximately the case in the regions 
where electromagnetic surface waves dominate the heat 
transfer. In the rest of the region, where near-field ra- 
diative transfer is non-rcsonant, increasing radius leads 
to increased contribution from propagating waves. The 
spectral conductance as the radius of the spheres is in- 
creased to 20 /im is shown in Fig. [8l The increased con- 
tribution from the non-resonant parts of the spectrum is 
evident from the graph. This has an important implica- 
tion from an experimental point of view. We should be 



careful not to increase the size of the sphere to such an 
extent that the resonant radiative transfer is swamped by 
the non-resonant radiative transfer. For a 20 /im sphere, 
we can see that the increase in conductance as the gap 
decreases is still predominantly due to electromagnetic 
surface waves as shown in Fig. [51 In Fig. [51 the in- 
crease in spectral conductance as the gap is decreased 
from 2000 nm is plotted. Compared to Fig. [51 the signa- 
ture of electromagnetic surface waves is clearer from this 
plot. 
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FIG. 6: Plot of spectral conductance between two silica 
spheres of 5 /im radii at gaps of 250 nm, 500 nm, and 1000 
nm from 0.04 eV to 0.16 eV. The gaps have been chosen so as 
to maintain the same ratio of gap to radius as for the curves 
shown in Fig. [S] 
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FIG. 7: Spectral conductance of spheres of radii 1 /im, 2 /im, 
and 5 /tm at gaps of 100 nm, 200 nm, and 500 nm respectively. 
The value of gap to radius for all the curves is 0.1. 



In Fig. [TO]the spectral conductance at 0.061 eV (corre- 
sponding to the first peak in Fig. [5]) between two spheres 
is plotted as a function of gap for different values of the 
radii. The exponent of a power law fit (of the form 
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FIG. 8: Plot of spectral conductance between two silica 
spheres of 20 fim radii at gaps of 200 nm, 400 nm, and 2000 
nm from 0.04 eV to 0.16 eV. Unlike the case of plane-to-plane 
near-field radiative heat transfer, where the contribution from 
surface polaritons dominate, the conductance between the two 
spheres has comparable contributions from the resonant and 
non-resonant regions. 



FIG. 10: Plot of spectral conductance between two spheres of 
equal radii at 0.061 eV as a function of gap for various radii. 
For each of the curves, the spectral conductance is computed 
for the same values of radius/gap. The markers on each curves 
corresponds to radius/gap values of 100, 50, 25, and 10. 



C. Total conductance 




Frequency (eV) 

FIG. 9: The increase in spectral conductance from 2000 nm 
gap in Fig. |5]is shown here for gaps of 200 nm and 400 nm. 



y = Ax^) to the data points in Fig. [TO]is -1.001, -0.984, 
-0.9024, -0.7781,-0.6119 for radius 1 fim, 4 /im, 10 /im, 20 
/im, and 40 /im respectively. As the radius increases from 
1 /im to 40 /im, the slope of the curve decreases, indicat- 
ing an increased contribution from propagating waves. 
The behavior at smaller radii can be predicted from the 
variation of conductance with gap between two planes 
and the proximity force approximation, as discussed in 
Sectionini For larger diameters, the proximity force type 
approximation is seen to fail. 



The frequency limits for the calculation of conductance 
is taken to be 0.041 eV to 0.164 eV. The contribution to 
conductance from the rest of the frequency spectrum at 
300 K is not significant for spheres of smaller radii. How- 
ever, it is seen from Fig. [5]that frequencies below 0.04 eV 
contribute to the non-resonant heat transfer. This will 
not affect the increase in radiative transfer as the gap 
is decreased because that increase comes predominantly 
from the regions supporting surface waves. The total 
conductance for spheres up to radius of 5 /im are plotted 
against gap in Fig. [TlJ The slope of -6 (approximate) for 
spheres of radius 20 nm and 40 nm, which can be approxi- 
mated as point dipoles when the gap between the spheres 
is much larger than the radius, is in correspondence with 
the results from the dipole approximation. However, the 
dipole theory predicts that the conductance should flat- 
ten and reach a finite value as the gap decreases to zero, 
i.e., a slope of zero. What happens in fact is that the 
near-field effects begin contributing as the gap decreases 
and the slope in fact decreases from -6 to approximately 
-1, once again corresponding to the asymptotic theory. 
However, as the radius increases to 5 /im, the slope de- 
creases further. This decrease is because of increased 
contribution from non-resonant regions of the spectrum, 
as seen in Fig. [S] For larger values of spheres, the con- 
ductance is plotted in Fig. [121 On a log-log scale as 
plotted in Fig. [T^ the near-field effects are apparent for 
the spheres of 1 /im and 2 /im. The reason it does not 
seem so for the spheres of larger diameter is because of 
the large contribution from the non-resonant parts of the 
spectrum, as evidenced from the curve corresponding to 
the results of classical radiative transfer for the 20 /im 
sphere. To understand the effects of near-field transfer 
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FIG. 11: Total conductance plotted as a function of gap. The 
number next to each curve is the corresponding radius of the 
spheres. For 20 nm and 40 nm spheres, the slope of the con- 
ductance vs gap curve is approximately -6 for values of gap 
larger than the radius of spheres. As the radius of the spheres 
is increased, the slope gradually changes to approximately -1, 
corresponding to the asymptotic theory. 



for larger spheres, the total conductance for 20 fiia and 
25 fim spheres is plotted in Fig. [131 The data points 
are fit with a curve G — Aix~^ + A2X + A3, where x 
is the gap, Aix~" is the near-field contribution, A2X is 
the contribution due to non-resonant parts of the spec- 
trum (as well as any changes from classical effects due 
to increase in view factor). We see that the value of the 
exponent n is 0.5574 for the 20 fim sphere and 0.5035 for 
the 25 fim sphere. If the asymptotic theory we valid at 
these values of the gap, one would have expected a value 
of n = 1. It is expected that as the gap decreases, the 
conductance will approach the form predicted in Eq. [1] 
by the asymptotic theory (not confirmed in this paper 
due to computational restrictions). 



VII. REGIME MAP FOR THE TWO-SPHERE 
PROBLEM 

Unlike the two half-space problem [U, i, i, i, i, H, 0, Hi , 
the two sphere problem, when other relavant length scales 
such as the skin depth are unimportant or \\/s\ ~ 1 (as in 
the case of silica spheres), where e is the dielectric func- 
tion of the sphere, has three length scales - the wave- 
length in consideration A, the radius of the spheres a, 
and the gap between the spheres, x. Depending on the 
ratio of the length scales, different approximate theories 
can be used to predict the conductances. These regions 
in which different theories are applicable can be repre- 
sented on a regime map with the two axes representing 
two non-dimensional length scales, as shown in Fig. [TH 
and discussed below. 

1. a 3> A and x > X: In this case classical radiative 



FIG. 12: Variation of total conductance with gap for various 
sphere radii. 
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FIG. 13: Variation of total conductance for spheres of radii 
20 fim and 25 /im with gap. The computation has been re- 
stricted to a minimum gap of 200 nm because of the stringent 
numerical requirements of convergence discussed earlier. No- 
tice that the plot is no longer on a log-log scale. 



transfer can be employed. However, while near- 
field effects may not be important, interference ef- 
fects can become important. When a « A, diffrac- 
tion effects prevent the usage of classical radiative 
transfer for even emission from a single sphere. 

2. a <C A and a <^ x X: When the dipole moment 
of the particles is the dominant contributor to the 
radiative transfer, the point dipole approximation 
can be used. However, it should be mentioned that 
even though the gaps are numerically small, this is 
not a near-field effect in the sense discussed in this 
paper. 

3. a; ^ a < A: In this case, the conductance is sees to 
vary linearly with 1 /x and is indicative of the valid- 
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FIG. 14: Regime map for the two sphere problem. Radius of 
spheres is a, the gap between them is x, and the wavelength 
of radiation is A. The numerical solution to the two sphere 
problem is valid everywhere, limited only by computational 
restrictions. 



flat plate. Most importantly, we expect n to be a num- 
ber between and 1. 

The work is supported by an ONR grant (Grant 
No. N00014-03-1-0835) through University of California, 
Berkeley. 

APPENDIX A: PROPERTIES OF VECTOR 
SPHERICAL WAVES 



£ (S, <f) X Vl;^^*{6, <f)).fd^ = (Al) 
£ (^m 4>) X V}l>{e, 0)) . rdn = 5ij,5„,, (A2) 
£ {y^n^ ^) X ^p?^* ■ = -SipS^, (A3) 



ity of the proximity approximation in this regime. 

The numerical solution to the two sphere problem as dis- 
cussed in this paper is valid in all regions of the map, 
limited only by computational demands. 



VIII. CONCLUSION 

Near-field radiative transfer between two spheres has 
been investigated for the very first time using fluctua- 
tional electrodynamics formalism. Numerical results for 
spheres of equal radii have been presented and analyzed. 
It is seen that the proximity force approximation type 
theory for near-field radiative transfer is valid only when 
the radius and gap satisfy the condition x ^ a < X. For 
larger spheres, the proximity force approximation, which 
is widely used in calculating forces, is not valid for the 
gaps analyzed in this paper. The purpose of solving the 
two sphere problem is to extend the theoretical and nu- 
merical formulations of near-field radiative transfer to 
configurations of objects which can be tested experimen- 
tally. The solution to the two-sphere problem gives us 
an estimate of the value of radiative conductance one 
can expect from such an experiment. For microspheres 
with radii of 25 /im, we expect a radiative conductance 
between two silica spheres around 4.5 nWK~^ at a gap 
of 200 nm. From Fig. [131 we see that we can expect a 
conductance of 5.6 nWK~^ at a gap of 100 nm. If the 
experimental configuration is not a two-sphere configu- 
ration but a sphere adjacent to a flat plate, the results 
of this chapter can be used as a guide to estimating the 
conductance. We expect trends to be similar - that is 
we expect the increase in conductance to be of the form 
Ax~^ , where x is the gap between the sphere and the 



kK [ MW(fcX)-M«*(fc,Odr' = 

<5ip<5™,a5(fc>z[^) {kaa)cl^'^*{kaa)') (A4) 

iV/„\^(fcX).iVW*(fc.r;)dr' = 

JVa 

6ip6rrLqa'^(^klazl^^* {kaa)Cl^\kaa)j (A5) 
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